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Abstract 

We employ a nonrelativistic version of QCD (NRQCD) to study 
heavy gg-bound states in the lowest approximation without fine struc- 
ture. We use gluon configurations on a 16 3 x 48-lattice at (3 = 6.2 from 
the UKQCD collaboration. For a bare quark mass near that of the 
b-quark (Ma = 1.6) we obtain the bound state masses for the S, P 
and both types of D-waves. We also detect signals for two types of 
hybrids (ggg-states) . The results are sufficiently accurate to confirm 
that the values of the D-w&ve mass from both D-waves coincide thus 
indicating that the cubical invariance of the lattice is restored to full 
rotational invariance to a good approximation. 

We also study S and P-wave masses for values of the range of bare 
quark mass Ma = 1.0, 1.3, 1.6 & 1.9 . The results confirm the idea 
that the 5/P-splitting is relatively insensitive to the value of the bare 
quark mass. 



1 Introduction 



Heavy quark bound states (J /if), T) are not only intrinsically interesting but 
because of the nearly non-relativistic motion of the constituent quarks, offer 
the prospect of numerical calculations considerably more tractable than the 
standard calculations for light quarks |lj, @, [|, ||. Several theoretical pro- 
posals have been made || and a number of numerical calculations carried 

out 0,1- 

For quark masses M ^> Aq C d, the quark velocity v <C 1, so that its 
dynamics are well modelled by an effective field theory, nonrelativistic QCD 
(NRQCD), in which the kinetic energy term is nonrelativistic [Q, [J. Clearly 
this approximation will only be valid for momenta very much less than some 
cut-off A = a -1 (a is the lattice spacing) of the order of the heavy quark mass 
M. Furthermore the effective action will contain correction terms which serve 
both to implement relativistic effects such as the hyperfine structure splitting 
and to compensate for the effects of the high momentum states which have 
been removed from the theory. These may be reorganised into a power series 
in the characteristic quark velocity determined by the dynamics. In practice 
this series is truncated at some finite order in the quark velocity and the 
couplings tuned so as to reproduce the low energy behaviour of full QCD. 

These systems serve as an important testing ground for lattice QCD since 
finite volume effects are expected to be much less important than for light 
quark states. Furthermore the NRQCD approach avoids the problem of 
fermion doubling and allows a rapid calculation of the quark propagator 
as an initial value problem (since it is determined by a lattice difference 
equation which is first order in time) || |j. Issues to be decided concern 
the appropriate size of lattice and the optimum values for the couplings and 
masses in the "non-relativistic" hamiltonian in order to best describe the 
bound state. 

In this letter we extend the numerical study of heavy quark bound states 
in QCD to a higher value of the coupling constant j3 = 6.2, using QCD gauge 
field configurations on 16 3 x 48-lattice from the UKQCD collaboration. We 
have employed the lowest order approximation in which we retain only the 
kinetic energy term in the effective action and studied correlation functions 
corresponding to the ground state (S'-wave) and angular excitations (P & D) 
of the gg-system. We have also looked for hybrid ggg-states and obtained a 
clear though still noisy signal. We have examined the S & P- waves for other 
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(lower) bare masses for the quarks and detected trends in the bound state 
energies and level splittings consistent with the results of Davis and Thacker 

§■ 

The important problem of fine structure of the energy levels we have not 
tackled in this paper partly because of restrictions in computational facilities. 
We intend to explore these and other questions in future work. 
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2 Propagators and Operators 

The evolution equation for the heavy quark propagator takes the form 

G (x, t + 1) = U$ (x, t) [{l - G (x,t)j + 5 (x) 8 (t) (1) 

where x labels the spatial position and G (t) = for t < and n is the order 
of the update as discussed by Thacker and Davies @. The modified update 
is necessary for stability at certain values of the bare quark mass. We have 
used n = 2 & 3, as appropriate. In fact there appeared to be little sensitivity 
to the value of n in the results when sufficient stability was present. The 
kinetic operator Hq is defined as 

#o = — -^AMr (2) 

The covariant finite differences A + , A - are given by their usual expressions 
(we have suppressed all colour and spin indices) 

A+G (x, t) = Ui (x, t) G (x + i,t) - G (x) (3) 

ArG (:r, t) = G (x, t) -Uj(x- i, t) G (x - i, t) (4) 
The correlation functions that we have considered are of the form 

g i (x,t) = (O i (x,t)Ol(0,0)) (5) 

where the index i runs over the different lattice operators corresponding to 
the S, P, D and hybrid states. For the S wave Os (x, t) with quantum 
numbers J PC = _+ 

O s {x,t) =xHx,t)iP(x,t) (6) 

The fields if) (x, t) and x { x i t) represent the quark and antiquark with propa- 
gators G (x, t) and G^ (x, t) respectively. In the current approximation these 
are spin independent. The higher angular momentum states involve the cor- 
relations of two displaced quarks, the P operator with J PC = l + ~ being 

P {x,t) = ^xHx,t)Aiif)(x,t) (7) 

= ^(xH^t)(^(^t))-(^x(x,t)) j iP(x,t)) (8) 
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The choice of A (with A now a forward derivative) ensures that the state 
is orthogonal to the S state and yields a correlation function that is entirely 
real. The continuum D state splits into two states under the lattice group, 
call them D\ and D 2 . These correspond to the operators 



Dl (x, t) = ^x 1 (x, t) {a* A,} ip (x, t) 



1 /~2 

D2 (x,t) = -x T (x,t) (A* 



A, )^(x,t) 



(9) 



(10) 



The indices i and j correspond to two (distinct) spatial directions and the 
curly braces in D\ indicate a symmetrised combination is to be taken. The 
hybrid states we investigate fit conveniently into the same computational 
framework. They correspond to taking commutators of the lattice covariant 
derivatives and therefore are only non-trivial in the presence of a gauge field 
configuration. We have 



o Hl {x,t) = l - x ] 0M)( 



Ai, Aj 



+ 



Ai, Aj 



ip(x,t) 



(11) 



H2 (x,t) = (x,t) 



Ai, A,- 



Ai, Aj 



i> (x, t) 



(12) 



These states possess quantum numbers J PC equal to 1~ + and 1 respec- 
tively. Because we have not included fine structure splitting in our hamilto- 
nian the 1 state does not mix with the gg-state with the same quantum 
numbers. An application of Wick's theorem then allows all the correlation 
functions to be built up out of sums of (gauge invariant) loops on the lattice 
constructed from the quark propagators and appropriate link matrices. 



3 Results 

Our results are obtained from ten uncorrelated configurations on a 16 3 x 48 
quenched lattice at f3 — 6.2 . They were used to generate the NRQCD 
correlation functions with bare quark masses Ma = 1.0, 1.3, 1.6 & 1.9. 
The configurations were obtained using the 64 i860-node Meiko Computing 
Surface at Edinburgh. The update consisted of a cycle of 1 three-subgroup 
Cabibbo-Marinari heat-bath sweep followed by 5 over-relaxed sweeps. 
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In order to extract the maximum information from each configuration we 
used three starting timeslices spaced by 16 units in the temporal direction 
and evolved the propagators from 4 3 evenly spaced spatial starting points 
within each of these slices. In addition we performed an average over spatial 
directions for the higher angular momentum states. It was found that each of 
these procedures contributed to a reduction in the raggedness of the results 
from a particular configuration. Indeed the ground state S'-wave correlation 
function computed from a single configuration was already rather smooth 
even out to 48 time slices. The higher waves of course required more data. 
A total of 6 quark propagators were needed to reconstruct all the correla- 
tors of interest. The errors were assessed by treating the data from each 
configuration as statistically independent. 

To extract estimates for the masses of the states we utilised a two expo- 
nential fit with four parameters to the zero spatial momentum correlator for 
times t < 24. 



The basis for the fit was a simple least squares procedure based on the sta- 
tistical errors estimated from the data. The fitting procedure was carried 
out for a sequence of initial time slices. The numbers quoted are those ob- 
tained from the best fit consistent within one standard deviation with those 
from succeeding initial time slices. We further checked that fits to a single 
exponential at large times gave statistically consistent results. 

For the case Ma = 1.6 we used all 10 gauge field configurations. If we use 
a value for the lattice spacing extracted from string tension measurements, 
aT 1 = 2.7 Gev, to convert masses to physical units this corresponds to M = 
4.32, roughly the b-quark mass. The procedure worked well particularly for 
S and P- waves. The fits are shown in Fig. 1. Clearly for these states we are 
in good quantitative control of both the statistical and systematic errors. As 
a further check we fitted the S'-wave data out to 48 timeslices and found a 
stable plateau in the ground state mass. Figs. 2a & 2b show the results for 
the two D states. Here we used the two exponential fit only for t < 20 since 
this yielded a more stable mass estimate. As we can see from Table 1. the 
two D-state masses are equal to within the errors which is an indication of 
the restoration of full rotational invariance at large distance. 

The data and fits for the hybrid states are plotted in Figs. 3a & 3b . 
Clearly our data is much noisier for these states and we attempted only a 




(13) 
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rough single exponential fit for an intermediate range of time-slices. This 
yielded estimates for the masses in the vicinity of the P-wave mass for both 
hybrids (see Table 1.). Although the results are very much provisional, it 
is nevertheless encouraging to see a signal at this stage. Hybrids are of 
particular interest because they can only exist as a result of propagation 
effects in the gauge fields and cannot be realised directly in a simple quark 
potential model. We hope in the future to achieve better results for these 
and the more conventional states. This will require both higher statistics 
and the use of improved operators with a stronger overlap on the physical 
states. 

On converting our lattice value for the S/ P-splitting to physical units we 
obtain AEgp = 0.35(2) Gev which is quite close to the experimental number 
of 0.4 for the splitting between triplet S and P-states of the T-system. It 
is also close to the implied physical value for the splitting found in ref M. 
This supports the physical ideas underlying the method. We have further 
strengthened the basis for the calculation by investigating the dependence of 
S/ P-splitting on the bare quark mass with results at Ma = 1.0, 1.3, 1.6 & 1.9 
summarised in Table 2. These confirm the approximate independence of this 
quantity on the bare quark mass. It should be noted however that these 
results which were obtained for reasons of computational economy from only 
3 configurations, required the use of a third order update method in order to 
encompass the lowest mass Ma = 1 . These new results for Ma = 1.6 using 
a third order update are consistent with the original results using a second 
order method at this value of the bare quark mass. 

Finally in Figs. 4a & 4b, we show effective mass plots for the S and 
P-waves at Ma = 1.6 and in Figs. 5a & 5b the same for Ma = 1.0 . In both 
cases the effective mass for the S'-wave shows a very convincing plateau. 
The P-wave at Ma = 1.6 is rather convincing also while the P-wave at 
Ma = 1.0 is much noisier. It is a matter of judgement whether these results 
can be taken as exhibiting the oscillations for the effective mass encountered 
by Davies and Thacker ||. It should be noted however that our computed 
errors are somewhat larger than those claimed in ref |J . The greater noise in 
the results at the lower mass is we suspect, due to the fact that with a third 
order update the limit of stability is roughly Ma = 1.0 @ . This potential 
instability underlines the need for using higher order updates. It is therefore 
encouraging that when there is stability (Ma = 1.3, 1.6 & 1.9) the results 
appear to be insensitive to the order of the update. This is an important 
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point that should be more thoroughly checked. 



4 Conclusions 

We have investigated the application of the NRQCD method in a numerical 
simulation of the lowest order kinetic energy hamiltonian, thus neglecting 
fine structure splitting for the purposes of this paper. We found that the 
S/ P splitting was roughly independent of a) the bare mass of the quark and 
b) the order of the time update used to compute the quark Green's functions. 
This encourages belief in the the physical meaning of the method. At the 
highest mass value we obtained results sufficiently accurate to confirm that 
the two versions of the .D-wave were degenerate. This was true not only of 
the .D-wave masses themselves but also of the additional parameters in the 
two exponential fit that mimic the higher mass states that obviously affect 
the form of the Green's function at small time separation. The S/ D-split is 
approximately twice that of the S/P-w&ve split. 

Hybrid states in the quark model have long been of interest [§]. It is 
encouraging therefore that were able to obtain a signal for hybrid ggg-states. 
Although the signal was noisy it suggests that they exist at a mass compa- 
rable to that of the D-wave. This implies an excitation energy above the 
ground state of ~ .8 Gev, a value rather less than but not grossly out of line 
with that suggested by excited static potential calculations applied to the T 



system fTD, 11 



Calculations to check the dependence of physical results on the bare quark 
mass suggests that this is weak as is the dependence on the order of the 
update. The absolute values of the energies of the states investigated do 
depend on the value of the bare quark mass. We note a trend of decreasing 
S-w&ve mass with increasing quark mass revealed in Table 2. In fact our 
results suggest that M^ 1 is approximately linear in the bare quark mass M 
as one might intuitively expect. This is consistent with the results of ref 
M . However we have not explored the effects of mass renormalization in 
this context, important though that issue is. 

We believe that we have obtained encouraging results for the NRQCD 
method and that it is possible to extract good results from only a few config- 
urations. We intend to extend our calculations to other sizes of lattice and 
other values of the coupling constant. This is particularly important since 
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the lattice value of ~ .13 for the 5/P-splitting suggests T radii of ~ 6 — 8 
lattice units at /3 = 6.2 . This is consistent with the wavefunction radius 
measured in refH and predicted in other calculations |T3 . In order to en- 



sure that finite size effects do not influence our results at this value of (3 it 
is important to repeat our calculations on a larger lattice. However because 
we see reasonable agreement with the results of ref|| obtained with (3 = 6.0 
on our size of lattice we do not expect there to be a strong effect. Further 
improvements in extracting data from QCD configurations can be envisaged 
such as better statistical treatment of the data and the use of better wave 
functions in the measurement of the more difficult states as well as higher 
statistics from more configurations. 
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Table Captions 

1. Masses for the states measured with Ma = 1.6 using the second order 
update. 

2. The Ms,M P & SyP-splitting at various masses using the third order 
update. 

Figure Captions 

Fig. 1 S and P wave correlation functions at Ma = 1.6 with solid line 
fit (two exponential). 

Fig. 2a D\ wave correlation functions at Ma = 1.6 with solid line fit 
(two exponential). 

Fig. 2b D 2 wave correlation functions at Ma = 1.6 with solid line fit 
(two exponential). 

Fig. 3a Hi hybrid correlation functions at Ma = 1.6 with solid line fit 
(single exponential). 

Fig. 3b if 2 hybrid correlation functions at Ma = 1.6 with solid line fit 
(single exponential). 

Fig. 4a Effective mass plot for S-wave at Ma — 1.6 . 
Fig. 4b Effective mass plot for P-wave at Ma — 1.6 . 
Fig. 5a Effective mass plot for S-wave at Ma — 1.0 . 
Fig. 5b Effective mass plot for P-wave at Ma — 1.0 . 
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Tables 



Table 1. 



State 


M a 


M x a 


S 


1.098(2) 


0.32(3) 


P 


1.230(10) 


0.29(4) 


Dl 


1.442(19) 


0.66(22) 


D2 


1.415(31) 


0.68(22) 


HI 


1.42(11) 




H2 


1.37(8) 





Table 2. 



Ma 


M s a 


M P a 


S/P-split 


1.0 


1.38(1) 


1.51(6) 


0.13(7) 


1.3 


1.20(1) 


1.32(2) 


0.12(3) 


1.6 


1.086(9) 


1.22(3) 


0.14(4) 


1.9 


1.003(8) 


1.14(2) 


0.14(2) 
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